home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libconv / TRY / ornl_dfir1d.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  4.5 KB  |  130 lines
  1.       subroutine ornl_dfir1d(    in_put, incinp, i0_inp, n_inp,
  2.      $                firfil, incfir, i0_fir, n_fir,
  3.      $                output, incout, i0_out, n_out,
  4.      $                alpha, beta)
  5.       double precision in_put(*), firfil(*), output(*)
  6.       double precision alpha, beta
  7.       integer incinp, i0_inp, n_inp
  8.       integer incfir, i0_fir, n_fir
  9.       integer incout, i0_out, n_out
  10. c-----------------------------------------------------------------------------
  11. c
  12. c   dconv1d  performs a 1D convolution in the time domain :
  13. c    O(j) = Beta*O(j) + Alpha * Sum[ I(i) * F(j-i) ]
  14. c
  15. c-----------------------------------------------------------------------------
  16. c
  17. c   PARAMETERS:
  18. c
  19. c    in_put:    Pointer to FIRST sample of sequence "in_put"
  20. c    incinp:    Increment between two successive values of "in_put"
  21. c    i0_inp: Index of the first element of "in_put" 
  22. c    n_inp:    Number of samples of "in_put"
  23. c
  24. c    firfil:    Pointer to FIRST sample of sequence "firfil"
  25. c    incfir:    Increment between two successive values of "firfil"
  26. c    i0_fir: Index of the first element of "firfil" 
  27. c    i0_fir:    Number of samples of "firfil"
  28. c
  29. c    output:    Pointer to FIRST sample of sequence "output"
  30. c    incout:    Increment between two successive values of "output"
  31. c    i0_out: Index of the first element of "output" 
  32. c    n_out:    Number of samples of "output"
  33. c
  34. c    beta:    Scaling factor for the Output Sequence on Entry 
  35. c    alpha:    Scaling factor for the convolution
  36. c-----------------------------------------------------------------------------
  37. c
  38. c   PLEASE NOTE:
  39. c    Please note that the array pointers must all point to the first 
  40. c    element of the array "i0_inp", "i0_fir" and "i0_out".
  41. c     If "in_put" for example is defined as
  42. c        dimension in_put(-25:45)
  43. c    Then "dconv1d" must be called with the following parameters
  44. c        call dconv1d( in_put(-25),1,-25,45, ... )
  45. c
  46. c-----------------------------------------------------------------------------
  47. c
  48. c   USAGE:
  49. c    This module compute the result of the convolution in the "output" range
  50. c    padding with zeroes when needed.
  51. c
  52. c    With some precautions it is possible that 
  53. c        the Output sequence   OVERWRITE   the Input one.
  54. c    Then (e.g. i0_out == i0_inp), it is necessary that   i0_fil >= 0
  55. c    Said in Digital Signal Processing (DSP) jargon: "firfil" must be CAUSAL
  56. c
  57. c    In theory, an input sequence of "n_inp" samples starting at time "i0_inp", 
  58. c    filtered by a sequence of "n_fir" samples starting at time "i0_fir",
  59. c    will result in a new signal of (n_inp + n_fir - 1) non zero samples starting at
  60. c    time (i0_inp + i0_fir). We just compute here the values that fall in that range
  61. c    and zero the rest.
  62. c    This may be interesting, for example when filtering a sequence of N samples, with a
  63. c    symmetric filter of 2m+1 samples. If one wants only to compute the central N resulting
  64. c    samples, the following call can be used:
  65. c        call dfir1d( f, 0, 1, N,  g, -m, 1, 2*m+1,  h, 0, 1, N)
  66. c
  67. c-----------------------------------------------------------------------------
  68. c
  69. c   This Fortran Subroutine written by
  70. c    Jean-Pierre Panziera
  71. c    Silicon Graphics Inc
  72. c    September 20, 1991
  73. c
  74. c-----------------------------------------------------------------------------
  75.  
  76.       double precision tmp, zero, one
  77.       parameter (zero = 0.0, one = 1.0)
  78.       integer i1_inp, i1_fir, i1_out, k0, k1
  79.       integer i, j, kout, kinp, kfir, n0, n1, jmin
  80. c############################################################################-
  81. c
  82. c
  83. c   Compute the FIR convolution
  84. c
  85. c
  86. c############################################################################-
  87.     i1_inp = i0_inp + n_inp - 1
  88.     i1_fir = i0_fir + n_fir - 1
  89.     i1_out = i0_out + n_out - 1
  90.     k0 = i0_inp + i0_fir
  91.     k1 = i1_inp + i1_fir
  92.  
  93.     n0 = i0_out
  94.     if( n0 .lt. k0)        n0 = k0
  95.     n1 = i1_out
  96.     if( n1 .gt. k1)        n1 = k1
  97. c-----------------------------------------------------------------------------
  98. c
  99. c    Scale the output
  100. c
  101. c-----------------------------------------------------------------------------
  102.     kout = 1 
  103.     do j = i0_out, i1_out
  104.         output(kout) = beta * output(kout)
  105.         kout = kout + incout
  106.     end do
  107. c-----------------------------------------------------------------------------
  108. c
  109. c    Compute the convolution
  110. c
  111. c-----------------------------------------------------------------------------
  112.     kout = 1 + (n1 - i0_out) * incout
  113.     do j = n1, n0, -1
  114.         tmp = zero
  115.         k0 = max(j- i1_fir, i0_inp)
  116.         k1 = min(j- i0_fir, i1_inp)
  117.         kinp = 1+ (k0-i0_inp) * incinp
  118.         kfir = 1+ (j-k0-i0_fir) * incfir
  119.         do i = k0, k1
  120.         tmp = tmp + in_put(kinp) * firfil(kfir)
  121.         kinp = kinp + incinp
  122.         kfir = kfir - incfir
  123.         end do
  124.         output(kout) = output(kout) + alpha * tmp
  125.         kout = kout - incout
  126.     end do
  127.  
  128.       return
  129.       end
  130.